############################################### File notes #######################################

############################### Clean workspace and define settings ##############################


### Run script to set up environment
rm(list=ls())


### Load packages if not already installed
install.packages("K:/BE/1261/Cases/Walgreens_Rite Aid_1610026/Research/Brian/Code/Code/R Packages/lfe_2.5-1998.zip", repos = NULL, type = "win.binary")
install.packages("K:/BE/1261/Cases/Walgreens_Rite Aid_1610026/Research/Brian/Code/R Packages/purrr_0.2.4.zip", repos = NULL, type = "win.binary")
install.packages("K:/BE/1261/Cases/Walgreens_Rite Aid_1610026/Research/Brian/Code/R Packages/tidyselect_0.2.3.zip", repos = NULL, type = "win.binary")
install.packages("K:/BE/1261/Cases/Walgreens_Rite Aid_1610026/Research/Brian/Code/R Packages/tidyr_0.7.2.zip", repos = NULL, type = "win.binary")
library(openxlsx)
library(foreign)
library(Matrix)
library(lfe)
library(data.table)
library(sp)
library(maps)
library(maptools)
library(ggplot2)
library(xtable)
library(lattice)
library(survival)
library(Formula)
library(Hmisc)
library(doBy)
library(readr)
library(tools)
library(tidyr)
library(readstata13)

dataIn <- file.path("K:","BE","1257","Projects","Medicare Claims","HospRetro","Drafts","Health_outcomes","FigureCreation","PSPatients","RefereeResponse","MakeTables")


#the filepath called "dataIn" had these .dta files in it (this is important for the code below):
# "death_fe_small.dta"          "death_lim_small.dta"         "diab_fe_small.dta"           "Diab_lim_small.dta"          "femAllBlocks.dta"            "fembetaBlockDiab.dta"        "fembetaDeath.dta"           
# "fembetaDiab.dta"             "fembetaDiab_overall.dta"     "fembetaHyp.dta"              "fembetaHyp_overall.dta"      "femblockDeath_overall.dta"   "femDiabBlocks.dta"           "femHypBlocks.dta"           
# "hyp_fe_small.dta"            "Hyp_lim.dta"                 "ll_death_small.dta"          "ll_diab_small.dta"           "ll_hyp_small.dta"            "maleAllBlocks.dta"           "malebetaBlockDeath.dta"     
# "malebetaBlockHyp.dta"        "malebetaDeath.dta"           "malebetaDiab.dta"            "malebetaDiab_overall.dta"    "malebetaHyp.dta"             "malebetaHyp_overall.dta"     "maleblockDeath_overall.dta" 
# "maleDiabBlocks.dta"          "maleHypBlocks.dta"           "mbm_death_small.dta"         "mbm_diab_small.dta"          "mbm_hyp_small.dta"           "overall_death_zip_small.dta" "Overall_diab_zip_small.dta" 
# "Overall_hyp_zip_small.dta"   "PS_YMeans_1_fillin.dta"      "PS_YMeans_2_fillin.dta"      "YMeans_fillin1.dta"          "YMeans_fillin2.dta"


############################################### create functions for later use ################ 
options(scipen = 999)
#Readin function
 fn.readin <- function(x){
   z <- read.dta13(paste0(dataIn,"/",x)) #dta. files
   z$dataname <- as.character(x) #Use the file name to determine sex results 
   z$sex <- as.character(x) #Use the file name to determine sex results 
   z$sex <- gsub(".dta","",gsub("_overall.dta","",z$sex)) 
   z$cond  <- gsub("fembeta","",gsub("malebeta","",z$sex)) #Use the file name to determine outcomes 
   z$cond  <- gsub("femblock","",gsub("maleblock","",z$cond)) #Use the file name to determine outcomes 
   z$sex  <- gsub("betahyp","",gsub("betadiab","",z$sex))
   z$sex  <- gsub("betaHyp","",gsub("betaDiab","",z$sex))
   z$sex  <- gsub("blockdeath","",gsub("betadeath","",z$sex))
   z$sex  <- gsub("blockDeath","",gsub("betaDeath","",z$sex))
   z$cond[z$cond=="hyp"]  <- "Hypertension"
   z$cond[z$cond=="Hyp"]  <- "Hypertension"
   z$cond[z$cond=="diab"]  <- "Diabetes"
   z$cond[z$cond=="Diab"]  <- "Diabetes"
   return(z)
 }

  
#BLOCKING FUNCTIONS 
#Limit dataset function
 fn.limit <- function(x){
   x <- x[,c("outcome","te","te_se","sex","cond")] #Keep relevant variables
   x <- x[!duplicated(x),] #Keep unique records
   x$te <- round(x$te,4) #Significant digits
   x$te[abs(x$te) > 2*x$te_se] <- paste0(x$te[abs(x$te) > 2*x$te_se] ,"*")
   x$te_se <- round(x$te_se,4)
   x$te_se <- paste0("(",x$te_se,")")
   
   return(x)
 }  
 #Limit to standard errors
 fn.se <- function(x){
   x <- x[,c("outcome","te_se","sex","cond")]
   names(x) <- c("outcome","te","sex","cond")
   x$order <-  1
   return(x)
 }  
 #limit to treatment effects
 fn.te <- function(x){
   x <- x[,c("outcome","te","sex","cond")]
   x$order <-  0
   return(x)
 }  
 #Final limitation
 fn.var <- function(x){
   x$Outcome[x$order==1] <- ""
   x <- x[,c("Outcome","ols.fem","ols.male","block.fem","block.male","me.fem","me.male")]
   names(x) <- c("Outcome","ols.Women","ols.Men","block.Women","block.Men","mat.Women","mat.Men")
   return(x)
 }  

#OLS FUNCTIONS
 #Limit ols 
 fn.limols <- function(x){
   names(x)[grepl("rndiabcomp1_",names(x))] <- "rn"  
   names(x)[grepl("rnami_",names(x))] <- "rn"  
   x <- x[x$rn == "1.postmerger",]
   x <- x[,names(x)[!grepl("rn",names(x)) & !grepl("dataname",names(x))]]
   name <- names(x)[grepl("sesq",names(x)) | grepl("beta",names(x))]
   x <- reshape(data = x, direction = "long", timevar="Outcome", idvar = c("sex","cond"), varying = names(x)[grepl("sesq",names(x)) | grepl("beta",names(x))], v.names="te")
   x$Outcome <- name
   row.names(x) <- NULL
   return(x)
 }  
 
 fn.fmtols <- function(z) {
   z$order[grepl("beta",z$Outcome)] <- 0
   z$order[grepl("sesq",z$Outcome)] <- 1
   z$Outcome[grepl("diabcomp1",z$Outcome)] <- "Asymptomatic"
   z$Outcome[grepl("diabcomp2",z$Outcome)] <- "Symptomatic"
   z$Outcome[grepl("glaucoma",z$Outcome)] <- "Glaucoma"
   z$Outcome[grepl("death",z$Outcome)] <- "Mortality"

   z$Outcome[grepl("acuteicd8",z$Outcome)] <- "Acute cardiac"
   z$Outcome[grepl("ami",z$Outcome)] <- "AMI"
   z$Outcome[grepl("ihd",z$Outcome)] <- "Ischemic heart"
   
   z$te[z$order==1] <-sqrt(z$te[z$order==1])
   z$te <- round(z$te,4)
   
   z$SE[z$order==0] <- z$te[z$order==1]
    z$te[z$order==0 & abs(z$te) > 2*z$SE] <- paste0(z$te[z$order==0 & abs(z$te) > 2*z$SE],"*")
    z$te[z$order==1] <- paste0("(",z$te[z$order==1],")")
    z$SE <- NULL  
 return(z)   
 }

 
 #MATCHING FUNCTIONS
 #Limit  
 fn.limmat <- function(x){
   names(x)[grepl("rndiabcomp1_",names(x))] <- "rn"  
   names(x)[grepl("rnami_",names(x))] <- "rn"  
   x <- x[grepl("postmerger",x$rn ),]
   x <- x[,names(x)[!grepl("rn",names(x)) & !grepl("dataname",names(x))]]
   name <- names(x)[grepl("sesq",names(x)) | grepl("beta",names(x)) | grepl("mem",names(x)) | grepl("mef",names(x)) ]
   x <- x[,name]
   return(x)
 }   
 
 fn.obcount <- function(x){
   for(i in names(x)[grepl("rn",names(x))] )  {
     x$flag[x[,c(i)]=="" ]  <- 1
   }
   for(i in names(x)[grepl("beta",names(x))] )  {
     x[is.na(x[,c(i)]),c(i)]  <- 0
   }
   x <- x[x$flag ==1,]
   for(i in names(x)[grepl("beta",names(x))] )  {
     x  <- x[!is.na(x[,c(i)]),]
     x[,c(i)]  <- max(x[,c(i)])
   }
   x <- x[ ,names(x)[grepl("beta",names(x)) ] ]
   x <- x[!duplicated(x) ,]
   z <- names(x)
   x <- reshape(data=x, direction = "long", timevar = "Outcome",v.names = "N", idvar = row.names(x), varying = list(names(x)) )
   x$Outcome <- z
   names(x) <- c("Outcome","N","id")
   
   
   x$Cond[grepl("_3",x$Outcome)] <- "Diabetes"
   x$Cond[grepl("_2",x$Outcome)] <- "Hypertension"
   x$Cond[grepl("_1",x$Outcome)] <- "All"
   
   x$Outcome[grepl("diabcomp1",x$Outcome)] <- "Asymptomatic"
   x$Outcome[grepl("diabcomp2",x$Outcome)] <- "Symptomatic"
   x$Outcome[grepl("glaucoma",x$Outcome)] <- "Glaucoma"
   
   x$Outcome[grepl("ihd",x$Outcome)] <- "Ischemic heart"
   x$Outcome[grepl("ami",x$Outcome)] <- "AMI"
   x$Outcome[grepl("acuteicd8",x$Outcome)] <- "Acute cardiac"
   
   x$Outcome[grepl("death",x$Outcome)] <- "Mortality"
   
   x$z1 <- x$Cond
   x$z2 <- x$Outcome
   x$Outcome[x$z2 == "Mortality"] <- x$z1[x$z2 == "Mortality"]
   x$Cond[x$z2 == "Mortality"] <- "Mortality"
   x <- x[,names(x)[names(x)!="z1" & names(x)!="z2"]]
   x <- x[!grepl("beta",x$Outcome),]
   x <- x[!(grepl("All",x$Cond) & grepl("AMI",x$Outcome)),]
   
   return(x)
   
 }   
 
 fn.rsq <- function(x){
   for(i in names(x)[grepl("rn",names(x))] )  {
     x$flag[x[,c(i)]=="rsq" ]  <- 1
   }
   x <- x[x$flag ==1 & !is.na(x$flag),]
   
   for(i in names(x)[grepl("rn",names(x))] )  {
     j <- paste0("beta_",gsub("rn","",i))
     x[,c(j)]  <- abs(x[,c(j)])
     x[,c(j)]  <- x[x[,c(i)]=="rsq" ,c(j)] 
   }
   
   
   x <- x[ ,names(x)[grepl("beta",names(x)) ] ]
   x <- x[!duplicated(x) ,]
   z <- names(x)
   x <- reshape(data=x, direction = "long", timevar = "Outcome",v.names = "N", idvar = row.names(x), varying = list(names(x)) )
   x$Outcome <- z
   names(x) <- c("Outcome","Rsq","id")
   
   
   x$Cond[grepl("_3",x$Outcome)] <- "Diabetes"
   x$Cond[grepl("_2",x$Outcome)] <- "Hypertension"
   x$Cond[grepl("_1",x$Outcome)] <- "All"
   
   x$Outcome[grepl("diabcomp1",x$Outcome)] <- "Asymptomatic"
   x$Outcome[grepl("diabcomp2",x$Outcome)] <- "Symptomatic"
   x$Outcome[grepl("glaucoma",x$Outcome)] <- "Glaucoma"
   
   x$Outcome[grepl("ihd",x$Outcome)] <- "Ischemic heart"
   x$Outcome[grepl("ami",x$Outcome)] <- "AMI"
   x$Outcome[grepl("acuteicd8",x$Outcome)] <- "Acute cardiac"
   
   x$Outcome[grepl("death",x$Outcome)] <- "Mortality"
   
   x$z1 <- x$Cond
   x$z2 <- x$Outcome
   x$Outcome[x$z2 == "Mortality"] <- x$z1[x$z2 == "Mortality"]
   x$Cond[x$z2 == "Mortality"] <- "Mortality"
   x <- x[,names(x)[names(x)!="z1" & names(x)!="z2"]]
   x <- x[!grepl("beta",x$Outcome),]
   x <- x[!(grepl("All",x$Cond) & grepl("AMI",x$Outcome)),]
   
   return(x)
   
 }   
 
 
 
 fn.mef <- function(x) {
   x <- x[,names(x)[grepl("mef",names(x))]]
   x <- x[!is.na(x[,1]),]
   name <- names(x)
   x$id <- row.names(x) 
   x <- reshape(data=x, direction = "long", timevar = "Outcome",v.names = "te.Women", idvar = "id", varying = list(names(x)[!grepl("id",names(x))]))
   x <- x[order(x$id,x$Outcome),]
   x$Outcome <- c(name,name)
   x$Outcome <- gsub("mef_","",x$Outcome)  
   x$Cond[grepl("_2", x$Outcome )] <-"Hypertension"  
   x$Cond[grepl("_3", x$Outcome )] <-"Diabetes"  
   x$te.Women[x$id==2] <- sqrt(x$te.Women[x$id==2])
   
   return(x)
 }
 
 fn.mem <- function(x) {
   x <- x[,names(x)[grepl("mem",names(x))]]
   x <- x[!is.na(x[,1]),]
   name <- names(x)
   x$id <- row.names(x) 
   x <- reshape(data=x, direction = "long", timevar = "Outcome",v.names = "te.Men", idvar = "id", varying = list(names(x)[!grepl("id",names(x))]))
   x <- x[order(x$id,x$Outcome),]
   x$Outcome <- c(name,name)
   x$Outcome <- gsub("mem_","",x$Outcome)  
   x$Cond[grepl("_2", x$Outcome )] <-"Hypertension"  
   x$Cond[grepl("_3", x$Outcome )] <-"Diabetes"  
   x$te.Men[x$id==2] <- sqrt(x$te.Men[x$id==2])
   return(x)
 }
 
 fn.beta <- function(x) {
   x <- x[,names(x)[grepl("beta",names(x))]]
   x <- x[x[,1] !=0,]
   name <- names(x)
   x$sex <- row.names(x) 
   x <- reshape(data=x, direction = "long", idvar = "sex",timevar = "Outcome",v.names = "te", varying = list(names(x)[!grepl("sex",names(x))]))
   x <- x[order(x$sex,x$Outcome),]
   x$Outcome <- c(name,name)
   x$sex[x$sex==3] <- "Men"
   x$sex[x$sex==4] <- "Women"
   x <- reshape(data=x, direction = "wide", timevar = "sex" , idvar = "Outcome")
   x$id <- "1"
   return(x)
 }
 
 fn.sesq <- function(x) {
   x <- x[,names(x)[grepl("sesq",names(x))]]
   x <- x[x[,1] !=0,]
   name <- names(x)
   x$sex <- row.names(x) 
   x <- reshape(data=x, direction = "long", timevar = "Outcome",v.names = "te", idvar = "sex", varying = list(names(x)[!grepl("sex",names(x))]))
   x <- x[order(x$sex,x$Outcome),]
   x$Outcome <- c(name,name)
   x$sex[x$sex==3] <- "Men"
   x$sex[x$sex==4] <- "Women"
   x$te <- sqrt(x$te)
   x <- reshape(data=x, direction = "wide", timevar = "sex" , idvar = "Outcome")
   x$id <- "2"
   
   return(x)
 }
 
 
 fn.matfmt <- function(x) {
   x$Cond[grepl("_3",x$Outcome)] <- "Diabetes"
   x$Cond[grepl("_2",x$Outcome)] <- "Hypertension"
   x$Cond[grepl("_1",x$Outcome)] <- "All"
   
   x$Outcome[grepl("diabcomp1",x$Outcome)] <- "Asymptomatic"
   x$Outcome[grepl("diabcomp2",x$Outcome)] <- "Symptomatic"
   x$Outcome[grepl("glaucoma",x$Outcome)] <- "Glaucoma"
   
   x$Outcome[grepl("ihd",x$Outcome)] <- "Ischemic heart"
   x$Outcome[grepl("ami",x$Outcome)] <- "AMI"
   x$Outcome[grepl("acuteicd8",x$Outcome)] <- "Acute cardiac"
   
   x$Outcome[grepl("death",x$Outcome)] <- "Mortality"
   
   z <- round(x[,names(x)[ grepl("Women",names(x)) | grepl("Men",names(x))]],4)
   x<- x[,names(x)[!grepl("te",names(x))]]
   x <- cbind(x ,z)
   x <- x[order(x$Cond,x$Outcome),]
   
   x <- x[ !grepl("beta",x$Outcome),]
   x <- x[ !grepl("sesq",x$Outcome),]
   x <- x[ !grepl("_",x$Outcome),]
   
   x <- x[!(x$Outcome=="AMI" & x$Cond=="All"),]
   
   x$z1 <- x$Cond
   x$z2 <- x$Outcome
   x$Outcome[x$z2 == "Mortality"] <- x$z1[x$z2 == "Mortality"]
   x$Cond[x$z2 == "Mortality"] <- "Mortality"
   x <- x[,names(x)[names(x)!="z1" & names(x)!="z2"]]
   x$z1 <- x$Cond
   x$z2 <- x$Outcome
   x <- x[,names(x)[names(x)!="z1" & names(x)!="z2"]]
   
   x <- x[order(x$Cond,x$Outcome,x$id),]
   x$m.se[x$id==1] <- x$te.Men[x$id==2]
   x$m.se[x$id==2] <- 0
   x$w.se[x$id==1] <- x$te.Women[x$id==2]
   x$w.se[x$id==2] <- 0
   
   x$te.Men[ x$id==1 & abs(x$te.Men)>=2*x$m.se ] <- paste0(x$te.Men[ x$id==1 & abs(x$te.Men)>=2*x$m.se ],"*")
   x$te.Women[ x$id==1 & abs(x$te.Women)>=2*x$w.se] <- paste0(x$te.Women[ x$id==1 & abs(x$te.Women)>=2*x$w.se],"*")
   
   x$te.Men[ x$id==2 ] <- paste0("(",x$te.Men[ x$id==2 ],")")
   x$te.Women[ x$id==2 ] <- paste0("(",x$te.Women[ x$id==2 ],")")
   x$ov[x$Cond=="Mortality"] <- 1
   x$ov[x$Cond=="Diabetes"] <- 2
   x$ov[x$Cond=="Hypertension"] <- 3
   x <- x[order(x$ov,x$Outcome),]
   x$m.se <- NULL
   x$w.se <- NULL
   return(x) 
 }

 
 ############################################### Execute functions and create tables ################ 
 
  
 #Execute the readin
 #Matching results
 #i keep the files labeled "death", only, because the file was created last and therefore has all of the results (ie., hypertension, death, diabetes) in it.
  files_match <- dir(dataIn)[(grepl("_lim",dir(dataIn)) | grepl("_zip",dir(dataIn)) | grepl("_fe",dir(dataIn))) & grepl("death",dir(dataIn))  ]
 
 match.data <- lapply(files_match , fn.readin )
 obscount <- lapply(match.data,fn.obcount)
 rsq <- lapply(match.data,fn.rsq)
 
 match.data <- lapply(match.data,fn.limmat) #limit the ols data 
 match.beta <- lapply(match.data,fn.beta)
 match.sesq <- lapply(match.data,fn.sesq)
 match.mef <- lapply(match.data[2:3],fn.mef)
 match.mem <- lapply(match.data[2:3],fn.mem)
 
 match.me <- lapply(1:length(match.mef), function(x) merge(match.mef[[x]],match.mem[[x]], by = intersect(names(match.mem[[x]]),names(match.mef[[x]]) ) , all= TRUE ))  
 
 match.coeff <- lapply(1:length(match.beta), function(x) rbind(match.beta[[x]],match.sesq[[x]]))  
 
 match.me <- lapply(match.me, fn.matfmt)
 match.me <- lapply(match.me,function(x) {names(x)[grepl("te",names(x))] <- gsub("te","me",names(x)[grepl("te",names(x))] )
 return(x) })
 
 match.coeff <- lapply(match.coeff, fn.matfmt)
 match.coeff <- lapply(match.coeff,function(x) {names(x)[grepl("te",names(x))] <- gsub("te","coeff",names(x)[grepl("te",names(x))] )
 return(x) })
 
 mat.me <- as.data.table(match.me[2])
 mat.me$id <- as.numeric(mat.me$id )-1
  names(mat.me)  <- c("order","Outcome","con","me.fem","me.male","ov") 
 
 #Diff-in-Diff results
  #ols results 
  
 files_ols <- dir(dataIn)[(grepl("fembeta",dir(dataIn)) | grepl("malebeta",dir(dataIn))) & !grepl("Block",dir(dataIn)) & !grepl("IHD",dir(dataIn)) & !grepl("ihd",dir(dataIn)) & !grepl("_overall",dir(dataIn)) & !grepl("block",dir(dataIn))]
 ols.data <- lapply(files_ols , fn.readin )
 ols.data <- lapply(ols.data,fn.limols) #limit the ols data 
 ols.data <- lapply(ols.data,fn.fmtols) #limit the ols data 
 
 olsdata <- rbindlist(ols.data) 
 olsdata <- reshape(olsdata, direction = "wide", idvar = c("cond","Outcome","order"), timevar = "sex")
 olsdata$con    <- olsdata$cond 
 olsdata$con[olsdata$Outcome == "Mortality"] <- "Mortality"
 olsdata$Outcome[olsdata$Outcome == "Mortality"] <- olsdata$cond[olsdata$Outcome == "Mortality"]
 olsdata <- olsdata[,c("Outcome","cond","order","con","te.fem","te.male")]
 names(olsdata) <- c("Outcome","cond","order","con","ols.fem","ols.male")
 olsdata$Outcome[olsdata$Outcome == "Death" ] <- "All"  
 
 #Blocking
 #blocking files identified by _overall
files_blocking <- dir(dataIn)[grepl("_overall",dir(dataIn))]
blocking <- lapply(files_blocking , fn.readin ) #readin the blocking data
blocking <- lapply(blocking,fn.limit) #limit the blocking data 
blockse <- lapply(blocking,fn.se) #SEs 
blockte <- lapply(blocking,fn.te) #TE 

block <- lapply(1:length(blocking),function(x) rbind(blockse[[x]],blockte[[x]])  )
block <- lapply(block,function(x) x[order(x$outcome,x$order),]  )
block <- rbindlist(block)
 block$outcome[block$outcome=="diabcomp1"]    <- "Asymptomatic" 
 block$outcome[block$outcome=="diabcomp2"]    <- "Symptomatic" 
 block$outcome[block$outcome=="ami"]    <- "AMI" 
 block$outcome[block$outcome=="ihd"]    <- "Ischemic heart" 
 block$outcome[block$outcome=="acuteicd8"]    <- "Acute cardiac" 
 block$outcome[block$outcome=="glaucoma"]    <- "Glaucoma" 
 block$outcome[block$outcome=="death"]    <- "Mortality" 
 block$con[block$outcome == "Mortality"] <- "Mortality"
 block$outcome[block$outcome == "Mortality"] <- block$cond[block$outcome == "Mortality"]
 names(block) <- c("Outcome","block","sex","cond","order","con")
 block <- reshape(block,direction = "wide", idvar = c("Outcome","order","con","cond"), timevar = "sex")
 block$Outcome[block$Outcome=="Death"] <- "All" 
 
ob <- merge(olsdata,block, by=intersect(names(olsdata),names(block)), all = TRUE)
ob <- merge(ob,mat.me, by=intersect(names(ob),names(mat.me)), all = TRUE)
ob <- ob[order(ob$cond,ob$ov,ob$Outcome,ob$order),] 
ob$con[ob$con=="Mortality"] <- ob$Outcome[ob$con=="Mortality"]
ob$Outcome[ob$con==ob$Outcome] <-  "Mortality"

obs <- split(ob,ob$con )
obs <- lapply(obs,  fn.var  )
obc <- rbindlist(obs)

addtorow <- list()
addtorow$pos <- list(0,0,1,2,10)
addtorow$command <- c( "& \\multicolumn{2}{c}{OLS} & \\multicolumn{2}{c}{Blocking}& \\multicolumn{2}{c}{Matching}\\\\\n", 
                       paste0(paste(gsub("mat.","",gsub("ols.","",gsub("block.","",names(obc)))),collapse = " & "),"\\\\\n"), 
                       paste0("& \\multicolumn{",length(obc)-1,"}{c}{\\underline{",names(obs[1]),"}} \\\\\n") , 
                       paste0("& \\multicolumn{",length(obc)-1,"}{c}{\\underline{",names(obs[2]),"}} \\\\\n") , 
                       paste0("& \\multicolumn{",length(obc)-1,"}{c}{\\underline{",names(obs[3]),"}} \\\\\n")) 
print(xtable(obc),add.to.row = addtorow,include.colnames = FALSE,include.rownames = FALSE, file=file.path(dataIn,paste0("Table_CompareEstimators.tex")))


